function out = curvelet_l1_min_BB( varargin )
% out = curvelet_l1_min_BB[Y,angles,truncate,imgSize,maxit,intialGuessStr,constThreshValue]
%
% Juergen Frikel 25/11/2011
%
% Sparse recovery: Accelerated iterated soft thresholding (AIST) with
%                  steplength selection and adapted universal threshold
%
% INPUT:
%
%   Y                   [noisy] projection Data 
%   
%   angles              measurement angles
%   
%   truncate            (optional) boolean flag. If its value is 'true', the dimensionality
%                       reduction in the curvelet domain is performed.
%
%   imgSize             (optional) size of the image to be reconstructed
%
%   maxit               (optional) max. number of iteration
%
%   initialGuessStr     (optional) initial guess
%
%   constThreshValue    regularization (thresholding) parameter
%
%
% OUTPUT variable 'out' is a struct array containing the following entries
%
%   rec                 reconstructed image
%
%   visibleCurveletIdx  indices of visible curvelets at the angular range
%                       given by 'angles'
%
%   iterates            cell array containing all iterates
%
%   error               vector of data errors of the corresponding iterates
%      
%
% Literature:
%   I. Loris, M. Bertero, C. De Mol, R. Zanella, L. Zanni. "Accelerating 
%   gradient projection methods for 1-constrained signal recovery by 
%   steplength selection rules", Appl. Comput. Harmon. Anal. 27 (2009),
%   247-254.


% fprintf(1,'Main Routine: Accelerated iterated soft thresholding (AIST)\n\n');
% fprintf(1,'\nProjection of CURVELET coefficients\n\n');


% intialize output variable
out = struct( 'rec',[],'visibleCurveletIdx',[],'iterates',[],'error',[],'rectime',[] );

% Main computation
tic;

[coeffs,out.visibleCurveletIdx,out.error,out.iterates,N] = gradien_descent_BB( varargin{:} );
% [coeffs,out.visibleCurveletIdx,out.error,out.iterates,N] = cg_iterate( varargin{:} );

out.rec     = ifdct( coeffs, out.visibleCurveletIdx, N, N);

totaltime = toc;
sec = round(totaltime);
s   = mod(sec,60) + round(1000*(totaltime-sec))/1000;
m   = (sec-mod(sec,60))/60;
out.rectime = sprintf('%2d min %2.3f sec\n\n',m,s);

%==========================================================================
function z = addCurveletCoeffs( x, factor, y, ccindex )
% compute: z = x + factor*y

lx = length(x);
ly = length(y);

if( lx~=ly )
    error('Coeff lengths must be equal!');
end

% allocate memory for the resulting coeff array
z = x;

% add approximation coeffs
if( size(x{1}{1})~=size(y{1}{1}) )
    error('Coeffs must have equal coeff sizes');
else
    z{1}{1} = x{1}{1} + factor.*y{1}{1};
end

for k=2:lx
   llx = length( x{k} );
   lly = length( y{k} );
   
   if( llx~=lly )
       error('Coeffs must have equal sizes');
   end
   
   for n=ccindex{k-1}
      if( size(x{k}{n})~=size(y{k}{n}) )
          error('Coeffs must have equal coeff sizes');
      end
      
      z{k}{n} = x{k}{n} + factor.*y{k}{n};
   end
end;
% z{end}{1} = x{end}{1} + factor.*y{end}{1};

%==========================================================================
function x = constMultiple( factor, x, ccindex )
% compute: z = x + factor*y

x{1}{1} = factor.*x{1}{1};

for k=2:length(x)
   for n=ccindex{k-1}
      x{k}{n} = factor.*x{k}{n};
   end
end;

%==========================================================================
function vector = extractCoeffs( coeffs, ccindex )
% extract visible curvelet coeffs to vector

% approximation coeffs
vector = coeffs{1}{1}(:);

% visible detail coeffs
for k=2:length(coeffs)
    for n=ccindex{k-1}
        vector = [vector;coeffs{k}{n}(:)];
    end
end

%vector = [vector;coeffs{end}{1}(:)];

%==========================================================================
function index = findVisibleCurveltIndices( maxscale, theta_min, theta_max )
% Input: 
%   maxscale    the finest (largest) scale for detail coeffs
%   theta_min
%   theta_max   angular range for the coeffs to be kept
%
% Output:
%   index       index array of DETAIL curvelet coeffs that should be kept
%               organization of the index array: 
%               index{1}    indices for coeffs{2}
%               index{2}    indices for coeffs{3}
%               ...
%               index{end}  indices for coeffs{end-1}

if( theta_max >= 180 )
    theta_max = 179;
end
theta_min = mod( theta_min, 180 );
theta_max = mod( theta_max, 180 );

theta_min = pi*theta_min/180;
theta_max = pi*theta_max/180;

if( maxscale < 3 )
    error('maxscale should be greater or equal than 3');
end

index = [];

tan_min = tan( theta_min );
tan_max = tan( theta_max );

if( tan_min>1e+10 )
    tan_min = Inf;
end
if( tan_max>1e+10 )
    tan_max = Inf;
end

itan_min = 1./tan_min;
itan_max = 1./tan_max;

for scale=3:maxscale

    N = 2*2^floor(scale/2);

    %fprintf(1,'tan_min = %f\n', tan_min);
    %fprintf(1,'tan_max = %f\n', tan_max);

    idx = [];
    
    idx_min = 0;
    idx_max = 0;

    dSlope = 2/N;
    % find idx_min
    if( tan_min <=1 && tan_min > -1 )
        for k=1:N
            slope1 = 1 - (k-1)*dSlope;
            slope2 = 1 - k*dSlope;

            if( tan_min<=slope1 && tan_min>slope2 )
                idx_min = k + N;
                break;
            end
        end
    else
        for k=1:N
           slope1 = -1 + (k-1)*dSlope;
           slope2 = -1 + k*dSlope;

           if( itan_min>=slope1 && itan_min<slope2 )
               idx_min = k;
               break;
           end
        end
    end
    % find idx_max
    if( tan_max <=1 && tan_max > -1 )
        for k=1:N
            slope1 = 1 - (k-1)*dSlope;
            slope2 = 1 - k*dSlope;

            if( tan_max<=slope1 && tan_max>slope2 )
                idx_max = k + N;
                break;
            end
        end
    else
        for k=1:N
           slope1 = -1 + (k-1)*dSlope;
           slope2 = -1 + k*dSlope;

           if( itan_max>=slope1 && itan_max<slope2 )
               idx_max = k;
               break;
           end
        end
    end

    % since the numbering of the curvelet coeffs starts at 135 degree, there
    % are 2 distinct regions of angles which have to be treated separately,
    % nameley [0 135] and [135 180] degrees.

    % chek if theta_min and theta_max are both within the same region
    border = 0.75*pi;
    sameregion = ( theta_min<=border && theta_max<=border ) ||...
                 ( theta_min>=border && theta_max>=border );

    tmp = idx_min;
    idx_min = min( [idx_min idx_max] );
    idx_max = max( [tmp idx_max] );
    
    % find indices of invisible curvelet-coefficients
    % Note: theta_min and theta_max describe the range of visible
    % curvelet-coeffs. Therefor, we are interested in the complement of the
    % index-set described by these angles
    if( theta_min<=theta_max)
        if( sameregion )
            % visible range is: [idx_min:idx_max]
            idx = idx_min:idx_max;
        else
            % visible range is: [1:idx_min idx_max:end]
            idx = [1:idx_min idx_max:2*N];
        end
    else
        if( sameregion )
            % visible range is: [1:idx_min idx_max:end]
            idx = [1:idx_min idx_max:2*N];
        else
            % visible range is: [idx_min:idx_max]
            idx = idx_min:idx_max;
        end
    end
    idx = mod( idx-1, 2*N ) + 1;

    
    % fill by symmetry
    index{scale-2} = [ idx idx+2*N ]; % index fuer Curvelet-Koeffs
end

index{end+1} = index{end};
%==========================================================================

function n = nofCurveletCoeffs( C, ccindex )

n = 0;
for k=1:length(C)
    if( k==1 )
        n = n + numel(C{1}{1});
    else
        for l=ccindex{k-1};%1:length(C{k})
            n = n + numel(C{k}{l});
        end
    end
end

%==========================================================================

function deg = degreeOfDimReduction( M,N,ccindex )

C = fdct_wrapping( ones(M,N) );
n = 0;
N = 0;
for k=1:length(C)
    if( k==1 )
        n = n + numel(C{1}{1});
        N = n;
    else
        for l=1:length(C{k})
            tmp = numel(C{k}{l});
            N = N + tmp;
            if( ismember(l,ccindex{k-1}) )
                n = n + tmp;
            end
        end
    end
end
deg = 1 - n/N;
deg = round(100*deg)/100;
% =========================================================================

function [wNew,ccindex,err,iterates,N] = gradien_descent_BB( varargin )

global exit_flag

% parse inputs
[Y,angles,truncate,imgSize,maxit,intialGuessStr,priorWeight] = parseinputs( varargin{:} );

% global Parameters for iteration
tol   = 0.01;   % tolerance

exit_flag    = false; % exit main routine on user 
plotIterates = true; % plot each iterate?
plotError    = false;
printLog     = true;
genericThreshold = false;
constThreshold = true;

% store all iterates in a 3D-array
iterates = [];%zeros(M,N,maxit);

N = imgSize;
if( isempty(N) )
    initialGuess = iradon(Y,angles);
    N = max(size(initialGuess));
end


switch lower(intialGuessStr)
    case 'zero'
        initialGuess = zeros(N);
    case 'fbp'
        initialGuess = iradon(Y,angles,N);
    case 'fbp_artefact_reduced'
        initialGuess = iradon_artefact_reduced(Y,angles,N);
    case 'backprojection'
        initialGuess = backproject(Y,angles,N);
    otherwise
        error('Unknown option for the initial guess!');
end

nofPr = size(Y,1);

% tv_weight = 0.001*(max(initialGuess(:))-min(initialGuess(:)));

% compute indices of visible coeffs
ccindex = [];
nbscales = ceil(log2(N) - 3);
if( truncate )
    %theta_min = 0;
    %theta_max = 90;
    
    % find curvlet indices corresponding to the angular range
    % [theta_min, theta_max]
    ccindex = findVisibleCurveltIndices( nbscales , min(angles(:)), max(angles(:)) );
    
    if( printLog )
        fprintf(1,'TRUNCATION of invisible curvelet coefficients\n');
        fprintf(1,'Visible angular range: [ %d, %d ]\n',min(angles(:)), max(angles(:)));
        fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
        fprintf(1,'MaxIterations = %d\n',maxit);
        fprintf(1,'Initial guess = %s\n',intialGuessStr);
    end
else
    % determine the full index array of all curvelet coefficients
    ccindex = findVisibleCurveltIndices( nbscales, 0, 179 );
    
    if( printLog )
        fprintf(1,'Computation of all coefficients\n'); 
        fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
    end
end

% Initializations
wOld = fdct( initialGuess, ccindex );
wNew = wOld;

dimRed = degreeOfDimReduction(N,N,ccindex);

threshSeq = [];
if( ~genericThreshold && ~constThreshold)
    threshSeq = getThresholdingSequence( wNew, ccindex, dimRed );
end
% if( genericThreshold )
%     [wNew,l1norm] = generic_softThresh_curvelet( wNew, ccindex, 0 );
% elseif( constThreshold )
%     [wNew,l1norm] = softThresh_curvelet_const( wNew, ccindex, constThreshValue );
% else
%     [wNew,l1norm] = softThresh_curvelet( wNew, ccindex, threshSeq );
% end


if( printLog )
    fprintf(1,'Dimension in the curvelet domain: %d\n',nofCurveletCoeffs(wOld,ccindex));
    fprintf(1,'Degree of DimReduction in the curvelet domain: %d\n',dimRed);
end

%%%%%%%%%%%%%%%%%%%%%%%%%
p = 1.01;
%%%%%%%%%%%%%%%%%%%%%%%%%
alpha   = 1000;        % initial inverse stepsize
E       = 10;       % parameter for the calculation of the acceptance cond.
gamma   = 1.e-4;    % parameter in the acceptance cond.
sigma1  = 0.1;      % parameters for finding new stepsize if the acceptance
sigma2  = 0.5;      % condition is 
epsilon = 1.e-10;   % 

% Barzilai-Borwein minimization
fprintf(1,'Main Routine: Barzilai and Borwein Minimization\n\n');
f    = zeros(E,1);
fval = minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
f(1) = fval;

[gOld, ng] = grad_minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
gNew = gOld;

% vector containing the data error of each iteration
err = [];

res = tol + 1;
it = 0;
if( plotIterates )
    figHandle = figure('KeyPressFcn',@keylistener);
    set(figHandle,'Position',[50 200 500 500])
    img = real( ifdct( wOld, ccindex, N, N ) );
    imagesc( img );
    colormap gray;
    title(strcat('Iteration ',num2str(it)));
    drawnow;
    
    if( plotError )
        plotHandle = figure;
        set(plotHandle,'Position',[600 200 500 500])
    end
end
while( res>=tol && it<maxit )
    it = it + 1;
   
   % check if the inverse stepsize 'alpha' is too small/large
   if( alpha <= epsilon || alpha >= 1/epsilon )
       alpha = setDelta( ng );
   end
   
   lambda = 1/alpha;
   
   wNew = addCurveletCoeffs( wOld, -lambda, gOld, ccindex );
   tmp_gOld = extractCoeffs( gOld, ccindex );
   g2   = tmp_gOld'*tmp_gOld;
   clear('tmp_gOld');
   
   tmp = max(f) + gamma*lambda*g2;
   fval  = minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
   
   % find stepsize that satisfy the acceptance condition
   zaehler = 0;
   while( fval > tmp && zaehler <50)
       sigma    = sigma1 + 0.5*(sigma2-sigma1);
       lambda   = sigma*lambda;
       
       wNew     = addCurveletCoeffs( wOld, -lambda, gOld, ccindex );%wOld - lambda*gOld;
       fval     = minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
       
       zaehler = zaehler + 1;
   end
   
   f(mod(it,E)+1) = fval;
   
   [gNew, ng] = grad_minFunc(Y,angles,nofPr,wNew,ccindex,p,N,priorWeight);
  
   if(it==1)
       dg = gOld;
   else
       dg = addCurveletCoeffs( gNew, -1, gOld, ccindex );
   end
   
   % Barzila-Borweit inverse stepsize
   tmp_gNew = extractCoeffs( gOld, ccindex );
   tmp_dg = extractCoeffs( dg, ccindex );
   
   alpha = -(tmp_gNew'*tmp_dg)/(lambda*g2);
   
   clear('tmp_gNew','tmp_dg');
   
   img = ifdct( wNew, ccindex, N, N );
   proj = radon( img, angles, nofPr );
   res  = norm( proj(:) - Y(:) )^2;
   err( it ) = sqrt(res);
   
   % positivity constraint
%    img(real(img)<0) = 0;
   
   wOld = fdct( img, ccindex);
   gOld = gNew;   
  
    if(mod(it,1)==0)
        if( plotIterates )
            %iterate = real( ifdct( wNew, ccindex, N, N ) );
            figure(figHandle);
%             img = real( ifdct( wNew, ccindex, N, N ) );
            imagesc( real(img) );
            colormap gray;
            title(strcat('Iteration ',num2str(it)));
            drawnow;
            %pause(1)
            
            if( plotError )
                figure(plotHandle)
                semilogy(it,res,'x');
                xlabel('iterations');
                ylabel('absolute error');
                hold on;
            end
        end
        if( printLog )
            fprintf(1,'%5d:  err=%e;  fval=%e;  norm(grad)=%e;  step=%e;\n',it,err(it),fval,ng,lambda);
        end
    end
%     
% %     if( mod(it,20)==0 )
% %         save( strcat('iterates',num2str(it),'.mat'),'iterates'); 
% %     end
%     
%     wOld = wNew;
%     %alpha = alpha_new;
%     %lambda_old = lambda;
%     
%     % save current iterate
%     iterates(:,:,it) = real( ifdct( wOld, ccindex, N, N ) )';
    
    if( exit_flag )
        return;
    end
end

%B = real( ifdct( wNew, ccindex, N, N ) );

if( printLog )
    %fprintf(1,'it=%5d;\terr=%e;\tfval=%e;\tl1-norm=%e;\tstep=%e;\tdiff=%e;\n',it,err(it),fval,l1norm,alpha,diff);
    fprintf(1,'Eof Main Computation\n\n');
end
%Eof

% =========================================================================


function [wNew,ccindex,err,iterates,N] = cg_iterate( varargin )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FUNKTIONIERT NOCH NICHT!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global exit_flag

% parse inputs
[Y,angles,truncate,imgSize,maxit,intialGuessStr,priorWeight] = parseinputs( varargin{:} );

% global Parameters for iteration
tol   = 0.01;   % tolerance

exit_flag    = false; % exit main routine on user 
plotIterates = true; % plot each iterate?
plotError    = false;
printLog     = true;
genericThreshold = false;
constThreshold = true;

% store all iterates in a 3D-array
iterates = [];%zeros(M,N,maxit);

N = imgSize;
if( isempty(N) )
    initialGuess = iradon(Y,angles);
    N = max(size(initialGuess));
end


switch lower(intialGuessStr)
    case 'zero'
        initialGuess = zeros(N);
    case 'fbp'
        initialGuess = iradon(Y,angles,N);
    case 'fbp_artefact_reduced'
        initialGuess = iradon_artefact_reduced(Y,angles,N);
    case 'backprojection'
        initialGuess = backproject(Y,angles,N);
    otherwise
        error('Unknown option for the initial guess!');
end

nofPr = size(Y,1);

% tv_weight = 0.001*(max(initialGuess(:))-min(initialGuess(:)));

% compute indices of visible coeffs
ccindex = [];
nbscales = ceil(log2(N) - 3);
if( truncate )
    %theta_min = 0;
    %theta_max = 90;
    
    % find curvlet indices corresponding to the angular range
    % [theta_min, theta_max]
    ccindex = findVisibleCurveltIndices( nbscales , min(angles(:)), max(angles(:)) );
    
    if( printLog )
        fprintf(1,'TRUNCATION of invisible curvelet coefficients\n');
        fprintf(1,'Visible angular range: [ %d, %d ]\n',min(angles(:)), max(angles(:)));
        fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
        fprintf(1,'MaxIterations = %d\n',maxit);
        fprintf(1,'Initial guess = %s\n',intialGuessStr);
    end
else
    % determine the full index array of all curvelet coefficients
    ccindex = findVisibleCurveltIndices( nbscales, 0, 179 );
    
    if( printLog )
        fprintf(1,'Computation of all coefficients\n'); 
        fprintf(1,'Nof DETAILS = %d\n',length(cell2mat(ccindex)));
    end
end

% Initializations
wOld = fdct( initialGuess, ccindex );
wNew = wOld;

dimRed = degreeOfDimReduction(N,N,ccindex);

threshSeq = [];
if( ~genericThreshold && ~constThreshold)
    threshSeq = getThresholdingSequence( wNew, ccindex, dimRed );
end
% if( genericThreshold )
%     [wNew,l1norm] = generic_softThresh_curvelet( wNew, ccindex, 0 );
% elseif( constThreshold )
%     [wNew,l1norm] = softThresh_curvelet_const( wNew, ccindex, constThreshValue );
% else
%     [wNew,l1norm] = softThresh_curvelet( wNew, ccindex, threshSeq );
% end


if( printLog )
    fprintf(1,'Dimension in the curvelet domain: %d\n',nofCurveletCoeffs(wOld,ccindex));
    fprintf(1,'Degree of DimReduction in the curvelet domain: %d\n',dimRed);
end

%%%%%%%%%%%%%%%%%%%%%%%%%
p = 1.05;
%%%%%%%%%%%%%%%%%%%%%%%%%
alpha   = 1;        % initial inverse stepsize
E       = 10;       % parameter for the calculation of the acceptance cond.
gamma   = 1.e-4;    % parameter in the acceptance cond.
sigma1  = 0.1;      % parameters for finding new stepsize if the acceptance
sigma2  = 0.5;      % condition is 
epsilon = 1.e-10;   % 

% Barzilai-Borwein minimization
fprintf(1,'Main Routine: Barzilai and Borwein Minimization\n\n');
f    = zeros(E,1);
fval = minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
f(1) = fval;

[gOld, ng] = grad_minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
gNew = gOld;
pOld = constMultiple( -1, gOld, ccindex );
pNew = pOld;

% vector containing the data error of each iteration
err = [];

res = tol + 1;
it = 0;
if( plotIterates )
    figHandle = figure('KeyPressFcn',@keylistener);
    set(figHandle,'Position',[50 200 500 500])
    img = real( ifdct( wOld, ccindex, N, N ) );
    imagesc( img );
    colormap gray;
    title(strcat('Iteration ',num2str(it)));
    drawnow;
    
    if( plotError )
        plotHandle = figure;
        set(plotHandle,'Position',[600 200 500 500])
    end
end
while( res>=tol && it<maxit )
    it = it + 1;
   
   % check if the inverse stepsize 'alpha' is too small/large
%    if( alpha <= epsilon || alpha >= 1/epsilon )
%        alpha = setDelta( ng );
%    end
   
   lambda = 1/alpha;
   
   wNew = addCurveletCoeffs( wOld, lambda, pOld, ccindex );
   tmp_pOld = extractCoeffs( pOld, ccindex );
   g2   = tmp_pOld'*tmp_pOld;
   clear('tmp_pOld');
   
   tmp = max(f) + gamma*lambda*g2;
   fval  = minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
   
   % find stepsize that satisfy the acceptance condition
%    zaehler = 0;
%    while( fval > tmp && zaehler <50)
%        sigma    = sigma1 + 0.5*(sigma2-sigma1);
%        lambda   = sigma*lambda;
%        
%        wNew     = addCurveletCoeffs( wOld, lambda, pOld, ccindex );%wOld - lambda*gOld;
%        fval     = minFunc( Y,angles,nofPr,wNew,ccindex,p,N,priorWeight );
%        
%        zaehler = zaehler + 1;
%    end
%    
   f(mod(it,E)+1) = fval;
   
   [gNew, ng] = grad_minFunc(Y,angles,nofPr,wNew,ccindex,p,N,priorWeight);
   
   % Fletcher-Reeves method
   tmp_gNew = extractCoeffs( gNew, ccindex );
   tmp_gOld = extractCoeffs( gOld, ccindex );
   
   n2_gNew = tmp_gNew'*tmp_gNew;
   n2_gOld = tmp_gOld'*tmp_gOld;
   beta = n2_gNew/n2_gOld;
%    clear('tmp_gOld','tmp_gNew');
   
   pNew = constMultiple( -1, gNew, ccindex );
   pNew = addCurveletCoeffs( pNew, beta, pOld, ccindex );
   
%    if(it==1)
%        dp = pOld;
%    else
       sigma = 0.01;
       wSig = addCurveletCoeffs( wOld, sigma, pOld, ccindex );
       [gSig, ~] = grad_minFunc(Y,angles,nofPr,wSig,ccindex,p,N,priorWeight);
       dg = addCurveletCoeffs( gSig, -1, gOld, ccindex );
%    end
   
   % Barzila-Borweit inverse stepsize
   tmp_dg = extractCoeffs( dg, ccindex );
%    tmp_pNew = extractCoeffs( pNew, ccindex );
   tmp_pOld = extractCoeffs( pOld, ccindex );
%    tmp_gNew = extractCoeffs( gNew, ccindex );
%    tmp_gOld = extractCoeffs( gOld, ccindex );
%    p2 = tmp_pNew'*tmp_pNew;
   
   alpha = (tmp_dg'*tmp_pOld)/(sigma*(tmp_gOld'*tmp_pOld));
%     alpha = -(tmp_pOld'*tmp_dp)/(beta*p2);

   clear('tmp_gOld','tmp_dg');
   
   img = ifdct( wNew, ccindex, N, N );
   proj = radon( img, angles, nofPr );
   res  = norm( proj(:) - Y(:) )^2;
   err( it ) = sqrt(res);
   
   % positivity constraint
%    img(real(img)<0) = 0;
   
   wOld = fdct( img, ccindex);
   gOld = gNew;   
   pOld = pNew;

    if(mod(it,1)==0)
        if( plotIterates )
            %iterate = real( ifdct( wNew, ccindex, N, N ) );
            figure(figHandle);
        %             img = real( ifdct( wNew, ccindex, N, N ) );
            imagesc( real(img) );%img(5:end-5,5:end-5)
            colormap gray;
            title(strcat('Iteration ',num2str(it)));
            drawnow;
            %pause(1)

            if( plotError )
                figure(plotHandle)
                semilogy(it,res,'x');
                xlabel('iterations');
                ylabel('absolute error');
                hold on;
            end
        end
        if( printLog )
            fprintf(1,'%5d: err=%e;  fval=%e;  norm(grad)=%e;  step=%e;  beta=%e;\n',it,err(it),fval,ng,lambda,beta);
        end
    end
%     
% %     if( mod(it,20)==0 )
% %         save( strcat('iterates',num2str(it),'.mat'),'iterates'); 
% %     end
%     
%     wOld = wNew;
%     %alpha = alpha_new;
%     %lambda_old = lambda;
%     
%     % save current iterate
%     iterates(:,:,it) = real( ifdct( wOld, ccindex, N, N ) )';
    
    if( exit_flag )
        return;
    end
end

%B = real( ifdct( wNew, ccindex, N, N ) );

if( printLog )
    %fprintf(1,'it=%5d;\terr=%e;\tfval=%e;\tl1-norm=%e;\tstep=%e;\tdiff=%e;\n',it,err(it),fval,l1norm,alpha,diff);
    fprintf(1,'Eof Main Computation\n\n');
end
%Eof

%% =========================================================================
function fval = minFunc(Y,angles,nofPr,coeffs,ccindex,p,N,priorWeight)
% This function implements the functional to be minimized

img = ifdct( coeffs, ccindex, N, N );
proj = radon(img,angles,nofPr);
res  = norm( proj(:) - Y(:) )^2;

% compute the weighted p-norm of the curvelet-coefficients
pnorm = 0;
for k=2:length(coeffs)
    for n=ccindex{k-1}
        pnorm = pnorm + sum( sum( abs(coeffs{k}{n}).^p ) );
    end
end

fval = res + priorWeight*pnorm^(1/p);

%% =========================================================================

function [grad, normgrad] = grad_minFunc(Y,angles,nofPr,coeffs,ccindex,p,N,priorWeight)
% This function implements the functional to be minimized

img = ifdct( coeffs, ccindex, N, N );
proj = radon( img, angles, nofPr );
proj = 2*backproject( proj - Y, angles, N );
grad = fdct( proj, ccindex );

% compute the weighted p-norm of the curvelet-coefficients
normgrad = 0;
pnormconst = (curvelet_pnorm_fine_scales( coeffs, ccindex, p ))^(1-p);
for k=2:length(coeffs)
    for n=ccindex{k-1}
        grad{k}{n} = grad{k}{n} + pnormconst.*priorWeight.*sign(coeffs{k}{n}).*abs(coeffs{k}{n}).^(p-1);
        normgrad = normgrad + sum( sum( abs(coeffs{k}{n}).^2 ) );
    end
end

normgrad = normgrad + sum( sum( abs(coeffs{1}{1}).^2 ) );
normgrad = sqrt( normgrad );

%% ==========================================================================

function pnorm = curvelet_pnorm_fine_scales( coeffs, ccindex, p )
pnorm = 0;
for k=2:length(coeffs)
    for n=ccindex{k-1}
        pnorm = pnorm + sum( sum( abs(coeffs{k}{n}).^p ) );
    end
end

%% ==========================================================================

function [Y,angles,truncate,imgSize,maxit,intialGuessStr,constThreshValue] = parseinputs( varargin )

narginchk(2,7);

if( nargin>=1 ) 
    Y = varargin{1};
end

if( nargin>=2 ) 
    angles = varargin{2};
end

truncate = false;
if( nargin>=3 ) 
    truncate = varargin{3};
end

imgSize = [];
if( nargin>=4 ) 
    imgSize = varargin{4};
end

maxit = 50;
if( nargin>=5 ) 
    maxit = varargin{5};
end
    
intialGuessStr = 'zero';
if( nargin>=6 )
    if( ~ischar(varargin{6}) )
        error('Initial guess parameter must be a string!');
    end
    intialGuessStr = varargin{6};
end

constThreshValue = 1;
if( nargin>=7 )
    constThreshValue = varargin{7};
end

%==========================================================================
function delta = setDelta(ng)

if( ng > 1 )
    delta = 1;
elseif( ng <= 1 && ng >= 1.e-5 )
    delta = 1/ng;
else
    delta = 1.e+5;
end

% =========================================================================
function keylistener(src,eventdata)

global exit_flag

if strcmp(eventdata.Key,'x')
    exit_flag = true;
    disp('Exit Main routine');
end